OPEN 3 ACCESS Freely available online 



•0-PLOS I o-^E 



Selection of Reliable Biomarkers from PGR Array 
Analyses Using Relative Distance Computational Model: 
Methodology and Proof-of-Concept Study 

Chunsheng Liu*^ Hongyan Xu^ Slew Hong Lam, Zhiyuan Gong* 

Department of Biological Sciences, National University of Singapore, Singapore, Singapore 



Abstract 

It is increasingly evident about the difficulty to monitor chemical exposure through biomarkers as almost all the biomarkers 
so far proposed are not specific for any individual chemical. In this proof-of-concept study, adult male zebrafish {Danio rerio) 
were exposed to 5 or 25 |xg/L 17p-estradiol (E2), 100 \ig/L lindane, 5 nM 2,3,7,8-tetrachlorodibenzo-p-dioxln (TCDD) or 
15 mg/L arsenic for 96 h, and the expression profiles of 59 genes involved In 7 pathways plus 2 well characterized blomarker 
genes, vtgl {vitellogenin!) and cyplal {cytochrome P450 lAl), were examined. Relative distance (RD) computational model 
was developed to screen favorable genes and generate appropriate gene sets for the differentiation of chemicals/ 
concentrations selected. Our results demonstrated that the known biomarker genes were not always good candidates for 
the differentiation of pair of chemicals/concentrations, and other genes had higher potentials In some cases. Furthermore, 
the differentiation of 5 chemicals/concentrations examined were attainable using expression data of various gene sets, and 
the best combination was the set consisting of 50 genes; however, as few as two genes (e.g. vtgl and hspaS [heat shock 
protein 5]) were sufficient to differentiate the five chemical/concentration groups In the present test. These observations 
suggest that multi-parameter arrays should be more reliable for blomonitoring of chemical exposure than traditional 
biomarkers, and the RD computational model provides an effective tool for the selection of parameters and generation of 
parameter sets. 
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Introduction 

Increasing attention has been drawn to the wide occurrence of 
natural and man-made chemicals in the aquatic environment. 
Many chemicals can be bioaccumulated in the aquatic organisms 
and magnified in the food chains, thus threatening human health. 
The Minamata disease is a typical case, where methylmercury 
(MeHg) poisoning occurred in human due to the ingestion of fish 
and shellfish contaminated by MeHg [1]. Such scenarios have 
promoted researchers to develop early-warning methods for 
monitoring contaminants in the aquatic system through both 
chemical monitoring and blomonitoring. 

As new pollutants in the environment are emerging rapidly, it 
becomes increasingly unfeasible to monitor all contaminants in the 
environment. Since the presence of a foreign chemical in a 
segment of the environment does not always indicate adverse 
biological effects [2], it is important to combine chemical 
monitoring with the blomonitoring for a reliable environmental 
risk assessment. An ideal approach is to examine biological 
responses that can reflect the contaminants in the exposed 
organisms [2]. Under this concept, various biomarkers from fish 
have been proposed and used for blomonitoring aquatic contam- 
inants. However, most of biomarkers proposed were not specific 



for individual chemicals. For example, biomarker for estrogen, vtgl 
mRNA could be induced not only by the native female hormone, 
1 yp-estradiol (E2), but also by many other compounds that can 
interact with estrogen receptors, including many xenobiotics, such 
as lindane [3]. The expression of cyplal was up-regulated by 
2,3,7,8-tetrachlorodibenzo-/)-dioxin (TCDD) as well as by other 
chemicals such as arsenic in mice [4]. 

It has been demonstrated that exposure to single chemicals 
generated unique gene expression signature in experimental 
animals [5-9]. Therefore, a multi-parameter quantitative real- 
time PGR (qRT-PCR) array could be developed as a useful tool to 
difierentiate a complicated set of chemical groups. However, in 
previous studies, the parameters (genes) were selected only based 
on responsive difference of gene expression among chemicals after 
exposure [10-1 1] and did not represent the best parameter (gene) 
set for the discrimination of chemicals. Therefore, a proof-of- 
concept study was designed and conducted in the present study, 
with the objective of finding the best parameter (gene) set for the 
discrimination of chemicals tested. Especially, a relative distance 
(RD) computational model was developed to select gene sets from 
61 gene examined for chemical discrimination. Therefore, it is 
feasible to integrate qRT-PCR arrays and RD computational 
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model to develop a reliable biomonitoring tool for chemical 
exposure. 

Materials and Methods 

Chemicals and reagents 

E2, lindane, TCDD and arsenic {Na2HAs04-7H20) were 
purchased from Sigma (St. Louis, MO, USA). Arsenic was 
dissolved in deionized water directly and the other three chemicals 
were dissolved in dimethyl sulfoxide (DMSO) as stock solutions. 
The TRizol reagent and LightCycle FastStart DNA Master SYBR 
Green I were obtained from Invitrogen (New Jersey, NJ, USA) and 
Roche Applied Science (Mannheim, Germany), respectively. 

Fish and chemical exposure 

In this study, experimental procedures were carried out 
following the approved protocol by Institutional Animal Care 
and Use Committee of National University of Singapore (Protocol 
079/07). Adult male zebrafish [Danio rerio, 5-month old) were 
purchased from a local aquarium farm (Mainland Tropical Fish 
Farm, Singapore), and acclimated for at least two weeks in our 
aquarium before chemical treatment. After acclimation, fish were 
exposed to 5 nM TCDD, 5 |xg/L E2, 50 |J,g/L E2, 100 |J,g/L 
hndane or 15 mg/L arsenic for 96 h in a static condition. Each 
tank (5 L size) included 3 L exposure solution and 3 fish, and each 
concentration included 3 replicated tanks. During the exposure 
period, fish were fed once a day with commercial frozen 
bloodworms (Ilikari) as described before [12]. The concentrations 
of these chemicals were chosen based on previous studies of ours 
and others [12-16], where biological effects ofthe.se concentra- 
tions have been confirmed by significant changes of some mRNAs 
examined. For E2, two concentrations were used to test the 
feasibilit)' to develop a gene expression based model to diflFeren- 
tiate exposure concentrations besides different chemicals. Fresh 
chemical solutions were daily replaced during the exposure 
experiment. For E2, lindane and TCDD exposure experiments, 
treatment and control groups received 0.01% DMSO, and for 
arsenic exposure experiments, treatment and control groups 
received 0.01% deionized water in this study. After 96-h exposure, 
the fish were anesthetized with MS-222 (1 mM) and livers were 
collected and preserved in TRizol reagent at -80°C until RNA 
isolation. 

Selection of target genes for PCR array 

A PCR array of sixty-one zebrafish genes was designed as 
follows. First, seven well characterized pathways commonly 
affected by chemicals were selected: oxidative and metabolic 
stress [17-18], apoptosis signaling [19-20], proliferation and 
carcinogenesis [21-22], DNA damage and repair [23-24], growth 
arrest and senescence [25-26], heat shock [27-28], and inflam- 
mation pathways [29-30]. Representative genes from these 
pathways were selected by referring Molecular Toxicology 
PathwayFinder PCR array from SABioscience Gene Network 
Central (http://www.sabiosciences.com/rt_pcr_product/HTML/ 
PAHS-3401Z.html). Second, annotated zebrafish orthologues of 
human genes were searched from Ensemble website and 
confirmed using online synteny tool [31]; unannotated zebrafish 
orhologues were manually determined first by amino acid 
sequence comparison with human candidate sequences through 
UCSC website (http://genome.ucsc.edu/) and then confirmed by 
comparison of genomic organization, chromosomal locations and 
chromosomal synteny analysis as conducted in a previously study 
[32]. Finally the zebrafish orthologues of 59 human genes were 
obtained for designing of PCR primers. In addition, two well- 



established biomarker genes, vtgl and cyplal, were also included in 
order to compare the potentials of biomonitoring between 
traditional biomarkers and genes/gene sets developed in this 

study, as inducers of vtgl and cyplal such as E2 and TCDD were 
also used in the present exposure experiments. The complete list of 
genes in PCR array and their PCR primeer sequences are 
presented in Table SI. The number of genes in each pathway was 
14, 10, 10, 6, 4, 13 and 2 for oxidative and metabolic stress, 
apoptosis signaling, DNA damage and repair, proliferation and 
carcinogenesis, growth arrest and senescence, heat shock and 
inflammation pathways, respectively. 

Quantitative real-time PCR (qRT-PCR) 

Total RNA was isolated from zebrafish livers with TRizol 
reagent and used for cDNA synthesis. Real time qPCR was 
performed using the LightCycler system (Roche Apphed Science, 
Mannheim, Germany) with LightCycler FastStart DNA Master 
SYBR Green I following manufacturer's instruction. The primer 
sequences were designed using Primer 3 software (http://frodo.wi. 
mit.edu/as). The amplicon efficiencies of primers were >90%. 
Three housekeeping genes, fi-actin [heta-actiri], fi-2m [beta- 2 -micro- 
globulin) and rpllSa [ribosomal protein L13a), were used as internal 
control and the geometric means the expression of the three 
housekeeping genes were used as the normalized factor by 2"^^^*^' 
method. Each group included three biological replicates and each 
replicate included a pool of three fish. 

Statistical analysis 

Gene expression values were logarithmically transformed (log2) 
before statistical analysis. The homogeneity and normality of data 
were examined using the Kolmogorov-Smirnov and Levene's test, 
respectively. Statically significant differences between treatment 
and corresponding control groups were evaluated by ANOVA 
based on a^-value <0.05. Average linkage (p < 0.05) was used to 
examine the cluster relationships of different treatment groups 
based on mRNA expression profiles. The statistical analyses were 
performed using Kyplot Demo 3.0 software (Tokyo, Japan). 

Relative distance (RD) computational model 

The differentiation of two chemical/ concentration groups not 
only depends on Euchdean distance between the two groups but 
also depends on the distance among individual replicates within 
each group. In this study, the RD computational model was 
developed to quantitatively describe the potential that three 
biological replicates from group A can be diflferentiated from the 
three replicates in group B based on mRNA expression profiles 
(fold change), and RD between one replicate from group A 
treatment and three replicates from group B {rdan) 

rdaib=mdaib—mdaa — 1/2 X SDaib — 1/2 X SD^^ (I) 



mdau = idam+dam + daib3)/^ (2) 



mdaa = {daUl + 4jl«3)/2 (3) 
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(4) 
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(5) 



SDalb = \l {{da\b\ - mdaibf + (da\bi - mdaibf + (dam - mda\bf)l{i - 1) 



SDaa = \J {(dalal - mdaaf + (4la3 - mdaaf)j (2-1) 



(7) 



where J is the total number of genes examined; a and b are gene 
expression values in the treatment groups A and B, respectively; 
m«4/i, is the mean Euclidean distance between one biological 
replicate from treatment group A (ay) and three replicates from 
treatment group B [bi, b2, bs); mdaa is the mean Eu[:hdean distance 
between one biological replicate from treatment group A (aj) and 
other two biological repUcates from the same group (a^, 05); SD^n 
is the standard deviation of EucUdean distance between one 
biological replicate from treatment group A («/) and three 
replicates from treatment group B [h], b2, i;;); 5D„„ is the standard 
deviation of Euclidean distance between one biological replicate 
from treatment group A treatment and other two biological 
replicates from the same group; dam, daib2 and dans are the 
Euclidean distance between one biological replicate from treat- 
ment group A (fl/) and three rephcates from treatment group B (bj, 
b2, hj); d„i„2 and are the Euchdean distance of biological 

responses between one biological replicate from treatment group A 
{di) and other two biological replicates from the same group (flp, 
as)- 

In this study, first, we calculated all the RD values between two 
chemical treatment groups using expression data of individual 
genes. When all six RD values wcr(; >() for each pair of chemicals, 
it was considered that the gene could be used to differentiate the 
two chemicals/concentrations. The cluster analyses (average 
linkage) were performed using commercial software (Kyplot Demo 
3.0, Tokyo, Japan) (/7-value <0.05) to confirm the feasibility of RD 
model in screening g("n("s for tlie differentiation of chemical/ 
concentration treatments. Second, the mean RD values were 
Ccdculated to quantitatively compare the potentials of individual 
genes in differentiating two chemicals/concentrations. Finally, a 
C-language computational program (see Program SI) was edited 
for selecting genes and generating gene sets that could be used to 
differentiate all of five chemical/ concentration treatments simul- 
taneously using the RD model developed in this study, and 
maximum mean RD of each gene sets with the same amount of 
genes and the corresponding components of genes were outputted. 



Results 

Broad changes of gene expression patterns in the seven 
selected pathways in response to chemical insults 

Adult male zebrafish were treated with 5 nM TCDD, 5 Hg/L 
E2, 50 Hg/L E2, 100 Hg/L hndane or 15 mg/L arsenic for 
96 hours and no mortalities were observed throughout the 
exposure experiment. As shown in Figure 1 and Table S2, 
exposure to different chemicals led to different gene expression 
profiles. TCDD exposure significantiy down-regulated the expres- 
sion of most selected genes involved in the oxidative and metabolic 
stress, apoptosis signaling, DNA damage and repair, proliferation 
and carcinogenesis, growth arrest and senescence, heat shock and 
inflammation pathways, while the expression of cyplal, h.spa5 and 
hsp70 (heat shock protein 70-kDd) was among the highest up- 
regulated. Treatment with arsenic significantly altered the 
expression of most selected genes in the seven pathways, such as 
up-regulation of expression ofptgsl {prostaglandin-endoperoxide synthase 
1), cyplal and hsp90aal [heat shock protm 90, alpha, class A member 1, 
tandem duplicate 1), and down-regulation of hlpl [Bcl-XL-like protein 
1), tnfr [tumor necrosis factor receptor) and I'tgl. A significant up- 
regulation in the expression of vtgl was observed upon exposure to 
5 or 50 |ig/L E2, clearly showing estrogenic activity. Similar to 
TCDD, exposure to E2 (5 or 50 Hg/L) significantiy down- 
regulated the expression of most selected genes included in the 
seven pathways investigated. In contrast, exposure to lindane up- 
regulated the expression of most selected genes in the seven 
pathways; with exception of only few down-regulated genes, 
notably cdknla (cyclin-dependent kinase inhibitor lA, transcript variant 1) 
in the growth arrest and senescence pathway and Jmo5 {flavin 
containing monooxygenase 5) in the oxidative and metabohc stress 
pathway. 

Correlation of RD and potential differentiation of 
chemical treatment pairs 

Using an RD computational model, we calculated all of RD 
values between two chemical/concentration treatment groups 
based on expression fold change of individual genes and the results 
are presented in Figure 2 (see details in Table S3) for all of the 10 
possible chemical/ concentration pairs. The ability of each of the 
6 1 genes to discriminate the chemical/ concentration pairs was 
tested by the software Kyplot Demo 3.0 program and the findings 
are presented in Figure 2. There was a good correlation of the RD 
and the ability to discriminate pair of chemicals/ concentrations. 
AH the genes with top and high RD values were found to be able 
to discriminate pair of chemicals/concentrations. For example, the 
two best known biomarker genes, vtgl and cyplal, were able to 
discriminate eight of the ten pairs: TCDD/arsenic, TCDD/ 
E2_high, E2_high/lindane, E2_high/arsenic, TCDD/E2_low, 
TCDD/lindane, lindane/E2_low, and arsenic/E2_low. However, 
for the lindane/arsenic pair, cyplal could not be used to 
discriminate them, while for the E2_low/E2_high concentration 
pair, both vtgl and cyplal failed to discriminate them. Interest- 
ingly, vtgl and cyplal were not always among the top of the list 
based on the calculated RD. There were also many other genes 
(even with better RD) that could be also used to differentiate the 
corresponding pair of chemicals. 

Selection of discriminating gene sets based on RD 
computational model 

While it is relatively easy to discriminate a pair of chemical 
treatment groups based on expression data from one or few genes, 
it is more challenging to discriminate multiple treatment groups. 
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Figure 1. Gene expression profiles included in seven selected pathways in male zebrafish livers after exposure to 100 |ig/L lindane, 
5 nM 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD), 5 |ig/L 1 7p-estradiol (E2), 25 |ig/L E2, or 15 mg/L arsenic for 96 h. There were 3 
biological replicates, and each replicate were pooled from 3 fish. Gene expressions were expressed as fold change relative to the corresponding 
control. The full names of genes can be found in Tables SI or S2. 
doi:1 0.1 371 /journal.pone.0083954.g001 



In the current dataset, no single gene can be used to discriminate 
all of the five chemical/concentration groups. Thus, it was 
necessary to select a gene set for discriminating the chemical/ 
concentration groups. Here, we further explored the RD model to 
select best gene sets for diflFerentiating all of the five chemical/ 
concentration groups. RDs were computed for all possible gene 



combinations from one to 6 1 genes and the highest mean distances 
for gene sets from 1 to 65 genes are presented in Figure 3. For 
example, the 2-gene set of the highest mean RD was vtgl and hspaS 
with a value of 10.57 (Figure 3 and Table S4) and the two genes 
can be used to discriminate the five chemical/concentration 
groups perfectiy (Figure 4A). In comparison, the gene pair of best 




Figure 2. Mean Relative Distances (RDs) between two chemicals/concentration groups. (A) TCDD vs. Arsenic; (B) TCDD vs. E2_high; (C) 
E2_low vs. E2_high; (D) E2_high vs. Lindane; (E) E2_high vs. Arsenic; (F) TCDD vs. E2_low; (G) TCDD vs. Lindane; (H) Lindane vs. Arsenic; (I) Lindane vs. 
E2_low; (J) Arsenic vs. E2_low. Black arrows indicate the positions of vtgl, and red arrows indicate the positions of cyplal; White boxes indicate the 
positions of genes that did not pass the model test and could not be used to discriminate the corresponding two chemicals/concentrations; Pink 
boxes indicate the positions of genes that passed the model test and could be used to discriminate the corresponding two chemicals/concentrations. 
TCDD: 5 nM 2,3,7,8-tetrachlorodibenzo-p-dioxin; lindane: 100 jj.g/L lindane; arsenic: 15 mg/L arsenic; E2_low: 5 ng/L 17(3 estradiol; E2_high: 50 |ig/L 
17p-estradiol. The information of RDs and the corresponding genes can be found in Table S3. 
doi:10.1371/journal.pone.0083954.g002 
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known biomarkers, vtgl and cyplal, has a vakie of 10.33 and they 
could not correctly discriminate all of the five groups, particularly 
the two concentration groups of E2 treatment (Figure 4B). AU 
other gene sets (3 or more genes) of the highest mean RD were also 
capable of differentiating all the five chemical/concentration 
groups correcdy (Figure 3). In general, there was an increase of 
mean RD with the number of genes in gene sets and the maximal 
mean RD (19.153) was observed in the set with 50 genes, where 
chemicals were also completely differentiated, including different 
concentrations (Figure 4C). 

Discussion 

The environment is continuously loaded with natural and man- 
made chemicals, and the effects of contaminant exposure to 
human health have been extensively documented [33-37]. In 
general, adverse effects of contaminants at population levels in 
wildlife and human tend to be delayed; when the effects finally 
become clear, the destructive processes may have been beyond the 
point where it can be reversed by available remedial actions [2] . 
Therefore, various biomonitoring methods have been developed 
in the past few decades for the purpose of early warning. However, 
most of these methods focused on one or several biological 
parameters (e.g., biomarkers vitellogenins and cytochrome P450 
enzymes lAl) [38-43]. To search for more biomarker genes to 
predict chemical contamination, it is common to use high 
throughput and large scale analyses such as DNA microarray 
and more recently RNA-seq platform [8,12,44]. However, the 
methodology for selecting biomarkers from thousands of genes 
could be a great challenge. Here we performed a proof-of-concept 
study by selecting a handful of biomarker genes to develop a 
practical assay with the aid of RD computational model. 




0 I I ■ I I I I 

0 10 20 JO iO ^0 60 70 

NunberofGeiKsiiGtae Seti 

Figure 3. Maximum mean RD of gene sets with different 
numbers of genes among 5 chemicals/concentrations. Blacl< 
arrow indicates the position of gene set (50 genes), where maximum RD 
was achieved. White box indicates the position of gene set (1 gene) that 
did not pass the model test and could not be used to differentiate the 
corresponding five chemicals/concentrations; Pink boxes indicate the 
positions of gene sets that passed the model test and could be used to 
differentiate the corresponding five chemicals/concentrations. The 
information about maximum mean RDs and the corresponding 
components of genes can be found in Table S4. 
doi:1 0.1 371/journal.pone.0083954.g003 



Here four chemicals including E2, lindane, TCDD and arsenic 
were tested. Both E2 and hndane exposures caused up-regulation 
of hepatic vtgl expression; similarly, treatment with TCDD or 
arsenic showed up-regulation of cyplal expression. These obser- 
vations are consistent with previous studies [3-4,45], suggesting 
the effectiveness of these chemical exposure experiments. In 
general, exposure to different chemicals resulted in different gene 
expression profiles in the seven biological pathways examined. For 
example, both of E2 and lindane induced vtgl expression, but E2 
down-regulated the expression of essentially all of the selected 
genes in the seven pathways while lindane up-regulated the 
expression of most of these genes. Similarly, TCDD down- 
regulated the expression of most of genes and arsenic up-regulated 
many of the genes, especially in two pathways, oxidative_and_- 
metabolic_stress and DNA_damage_and_repair, suggesting a 
molecular basis for their discrimination. 

In the current data set, we found that none of the 61 genes could 
be used to correctly discriminate all of the five chemical/ 
concentration groups; thus, it has to rely on multiple gene sets 
for successful discrimination, which should be the direction for 
future development of multiple gene signatures for discrimination 
of a multiple chemical groups, as previously proposed [8,12]. To 
systematically select the best discriminator genes, here we 
developed a computational model using RD to determine the 
prediction power of each gene or in combination with others. First, 
we demonstrated that there was a positive correlation between the 
RD values and the cUscrimination of different treatments groups 
(Fig. 2). In our data set, a minimum of two genes (e.g. vtgl and 
hspaS) could be used to successfully discriminate all of the five 
chemical/concentration groups. There is a general increase of 
mean RD values with the number of genes added to the gene set, 
which indicate the power of using more genes for discriminating 
more complicated data set. In our dataset, we also found that the 
50-gene set had the highest mean RD values, indicating that there 
is an optimal gene number used for the discrimination. From a 
practical viewpoint, the used of minimal number of genes wiU 
minimize workload and ease downstream data analysis. However, 
using more genes, especially those representing different molecular 
pathways, provides additional important biological information in 
molecular-marker based biomonitoring. 

In summary, the data of this study demonstrated chemicals that 
induced similar responses in biomarker (e.g., TCDD and arsenic, 
E2 and hndane) could cause different biological responses 
depending on the parameters examined, and the use of parameter 
sets consisting of different biological responses for biomonitoring 
should be more appropriate. Furthermore, the computational 
model based on RD may be useful to select appropriate gene sets 
to develop efficient biomarker-based biomonitoring. Considering 
the rapid, sensitive, convenient and high-throughput properties of 
PCR, a PCR array including multiple gene parameters should be 
a feasible tool to develop for biomonitoring of chemical exposure. 

Supporting information 

Table SI Sequences of primers for selected genes. 

(DOC) 

Table S2 mRNA expression profiles in the livers of 
zebrafish after chemical exposure. 

(DOC) 

Table S3 Mean relative distances (MRDs) of individual 
genes between chemicals. 

(DOC) 
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(A) v^l and hspaS 



(B) cx'plal and vr^i 




Figure 4. Clustering relationships among chemicals/concentrations using mRNA expression data of (A) cyplal and vtgl, (B) vtgl and 
hspaS and (C) 50 genes with the marximum RD. TCDD: 5 nM 2,3,7,8-tetrachlorodibenzo-p-dioxin; lindane: 100 |ig/L lindane; arsenic: 15 mg/L 
arsenic; E2_low: 5 j^g/L 17p-estradiol; E2_high: 50 ng/L 17p-estradiol. The full names of genes can be found in Table SI or S2. 
doi:1 0.1 371 /journal.pone.0083954.g004 



Table S4 Maximum mean relative distances (MMRDs) 
of gene sets with different amounts of genes among 3 
chemicals/ concentrations and the corresponding com- 
ponents of genes. 

(DOC) 

Program SI 

(ZIP) 
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